An important geoscience data structure is a voxel (voxset), which is a 3D cell-based model of the subsurface volume of part of the Earth. Voxel models are commonly used to represent some property of the subsurface in which the entire volume of interest is subdivided into cells, each with a value that represents a property of the cell volume. A voxel property can also be a vector quantity, in which case the property of a voxel cell will have a vector quantity defined by (vx, vy, vz).
We work with voxels using the Vox class and the Vox_display class. The Vox class can be used to open a Geosoft voxel (or vector voxel) or create a new Geosoft voxel, and the Vox_display class can create a visualization of a voxel that can be placed in a 3D view.
In [1]:
import geosoft.gxpy.gx as gx
import geosoft.gxpy.view as gxview
import geosoft.gxpy.group as gxgroup
import geosoft.gxpy.vox as gxvox
import geosoft.gxpy.vox_display as gxvoxd
import geosoft.gxpy.viewer as gxviewer
import geosoft.gxpy.utility as gxu
import geosoft.gxpy.map as gxmap
import geosoft.gxpy.surface as gxsurf
from IPython.display import Image
import numpy as np
gxc = gx.GXpy()
url = 'https://github.com/GeosoftInc/gxpy/raw/9.3.1/examples/data/'
vox_file = 'rjsmith_voxi_density.geosoft_voxel'
gxu.url_retrieve(url + vox_file)
gxu.url_retrieve(url + vox_file + '.xml')
vectorvox_file = 'mvi.geosoft_vectorvoxel'
gxu.url_retrieve(url + vectorvox_file)
gxu.url_retrieve(url + vectorvox_file + '.xml')
Out[1]:
In [2]:
# open the vox
with gxvox.Vox.open(vox_file) as vox:
# show some of the vox properties
print(vox.name)
print('dimensions (nx, ny, nz):', vox.nx, vox.ny, vox.nz)
print('coordinate system:', vox.coordinate_system)
print('surface extent:', vox.extent_xy)
print('volume extent:', vox.extent)
print('origin location:', (vox.origin_x, vox.origin_y, vox.origin_z))
print('unit of measure:', vox.unit_of_measure)
print('depth cell sizes in', vox.coordinate_system.unit_of_measure, ', from the bottom:')
for cell in vox.cells_z:
print('\t', cell)
The vertical dimension of a voxel (the 'z' dimension), can be referenced by elevation or by depth. If referenced by elevation, which is the default, the positive z direction is up and values can be assumed to represent elevations relative to the vertical reference datum of the coordinate system. Also, when referenced by elevation, the voxel origin is on the bottom slice of the voxel and z arrays of voxel cells sizes or locations start at the bottom and sequence up, in the positive direction.
You can also choose to work with a voxel by depth, either by setting the depth=True parameter when a new voxel is created, or an existing voxel is openned, or by changing the is_depth property to True.
For example, here we change the reference to depth and print depth-based properties. Note the z property differences from the example above.
In [3]:
# open the vox
with gxvox.Vox.open(vox_file, depth=True) as vox:
# show some of the vox properties
print(vox.name, 'properties by depth...')
print('volume extent:', vox.extent)
print('origin location:', (vox.origin_x, vox.origin_y, vox.origin_z))
print('depth cell sizes in', vox.coordinate_system.unit_of_measure, ', from the top:')
for cell in vox.cells_z:
print('\t', cell)
A Vox instance can be displayed withion a 3D view either as a solid with each cell assigned a color, or as a set of iso-surfaces. Vector voxels can be displayed as a solid, as vectors or as iso-surfaces.
Here we create a solid with starting minimum value of 0.005, though this can be changed dynamically in the Geosoft 3D viewer. The figure_map() method will render the solid into a group of a 3d view in a 3d map.
In [4]:
# open the vox
with gxvox.Vox.open(vox_file) as vox:
# create a display instance
with gxvoxd.VoxDisplay.solid(vox) as voxd:
# limit the display to densities above 0.005 g/cc
voxd.shell_limits = (0.005, None)
# Create a title and include the coordinate system.
title = "VOXI Density Model - RJ Smith Test Range" + '\n' + str(vox.coordinate_system)
# create a figure map for the vox. This creates a 3D view that is placed on a map and annotated
# to create a default presentation figure. By not specifying a name for the map a temporary map
# is created only for the life of this Python session.
vox_map_file = voxd.figure_map(title=title, legend_label='density\ng/cc').file_name
# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))
Out[4]:
In [5]:
# create an isosurface dataset, shell at 0.005
isosurface_file = gxsurf.SurfaceDataset.vox_surface(gxvox.Vox.open(vox_file), 0.005, overwrite=True).file_name
# create a temporary 3d view and place the isosurface dataset in the view
with gxview.View_3d.new() as v3d:
v3d_file = v3d.file_name
gxgroup.surface_group_from_file(v3d, isosurface_file)
# display here
Image(gxmap.Map.open(v3d_file).image_file(pix_width=800))
Out[5]:
Here we create three isosurfaces with specified colors.
In [6]:
# create an isosurface dataset, shells at 0.002, 0.01 and 0.05
isosurface_file = gxsurf.SurfaceDataset.vox_surface(gxvox.Vox.open(vox_file),
surfaces=(0.002, 0.01, 0.05),
color=(gxgroup.C_YELLOW, gxgroup.C_GREEN, gxgroup.C_RED),
overwrite=True).file_name
# create a temporary 3d view and place the isosurface dataset in the view
with gxview.View_3d.new() as v3d:
v3d_file = v3d.file_name
gxgroup.surface_group_from_file(v3d, isosurface_file)
# display here
Image(gxmap.Map.open(v3d_file).image_file(pix_width=800))
Out[6]:
A simple way to create a new voxel is to collect the voxel data into a 3D numpy array in which each dimension represents the cells in the X, Y and Z spatial directions. Data in an array is organized by X, then by Y, and lastly by Z, which results in data indexing being in the reverse order. For example, a voxel that has spatial dimensions (nx, ny, nz) of (100, 65, 20) will be stored in a 3d array with shape (20, 65, 100).
Note also that the Z direction is positive up, which means that a voxel origin, which is the center of the first cell (cell[0, 0, 0]) is the at the center of the deepest, lower-left cell in the voxel.
In this example we will start with a simple indexed data array in which the value of each cell is the Z index.
In [7]:
# create a 3D data array in which each z slice is set to the z index
data = np.empty((20, 65, 100), dtype=np.float64)
for i in range(data.shape[0]):
data[i, :, :] = i
# create a new voxel, and display
with gxvox.Vox.new(name='my vox', data=data, overwrite=True) as vox:
vox_map_file = gxvoxd.VoxDisplay.solid(vox).figure_map().file_name
# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))
Out[7]:
In the last example we created a voxel without true spatial dimensions. In this example we will scale the same voxel to have true spatial dimensions and be located at a known location on Earth. For this example we will assume that the voxels are 50 m square cells and the model is located on the NAD83 datum using projection UTM zone 19N.
In [8]:
# create spatially located voxel
with gxvox.Vox.new(name='my vox',
data=data,
overwrite=True,
origin=(450025, 6800025, -975), # the origin is the center of the bottom cell[0, 0, 0]
cell_size=50,
coordinate_system='NAD83 / UTM zone 19N') as vox:
# display in a map
vox_map_file = gxvoxd.VoxDisplay.solid(vox).figure_map(title=str(vox.coordinate_system)).file_name
# display the map here
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))
Out[8]:
Voxel cell sizes may vary along any dimension. For example, when defining a voxel model for geophysical modelling one might increase the vertical size of voxel as a function of depth. To do this you can provide the cell dimensions as an array of dimensions along any, or all axis. In this example we will increase the cell size as a function of depth, with each deeper cell being 10% larger. We will referencing depth rather than elevation as it makes the cell depth calculation simpler.
In [9]:
# create a z cell dimension increasing 10% as a function of depth (positive down)
nz = data.shape[0]
cell = 50
z_cell = [cell] * nz
for i in range(1, nz):
z_cell[i] = cell * (1. + i * 0.1)
# create voxel, and display
with gxvox.Vox.new(name='my vox',
data=data,
overwrite=True,
depth=True, # work using z as depth
origin=(450025, 6800025, 25), # the origin is the center of the top cell when working from depth
cell_size=(cell, cell, z_cell), # the z_cell size array defines variable z cell sizes
coordinate_system='NAD83 / UTM zone 19N') as vox:
vox_map_file = gxvoxd.VoxDisplay.solid(vox).figure_map(title=str(vox.coordinate_system)).file_name
# display here
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))
Out[9]:
Voxel cells may also store a vector value defined by a (vx, vy, vz) for each cell. For example, Geosoft MVI models define model magnetization as a vector quantity. Vector voxel files have extension .geosoft_vectorvoxel.
The is_vectorvox property of a Vox instance that holds a vector voxel will be True, and the Numpy array returned by the Vox.np() method will have 4 dimensions with the last dimension being 3 which contains the (vx, vy, vz) vector for each voxel cell.
Displaying a vector voxel as a solid will display the Vox as scalar colour based on vector amplitude, as follows.
In [10]:
# open the vector voxel
with gxvox.Vox.open(vectorvox_file) as vox:
# create a solid display
with gxvoxd.VoxDisplay.solid(vox) as voxd:
# limit the display to MVI susceptibility above 0.01 SI
voxd.shell_limits = (0.04, None)
# Create a title and include the coordinate system.
title = "VOXI MVI Susceptibility" + '\n' + str(vox.coordinate_system)
# create a figure map for the vox.
vox_map_file = voxd.figure_map(title=title, legend_label='MVI SI').file_name
# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))
Out[10]:
In [11]:
# open the vector voxel
with gxvox.Vox.open(vectorvox_file) as vox:
# create a vector display
with gxvoxd.VoxDisplay.vector(vox) as voxd:
# limit the display to MVI susceptibility above 0.04 SI
voxd.shell_limits = (0.04, None)
# set the vector specifications
voxd.vector_cone_specs = (3, None, None, 1000)
# Create a title and include the coordinate system.
title = "VOXI MVI Susceptibility"
# create a figure map for the vox.
vox_map_file = voxd.figure_map(title=title, legend_label='MVI SI').file_name
# create a PNG from the figure map and display in Jupyter
Image(gxmap.Map.open(vox_map_file).image_file(pix_width=800))
Out[11]: